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' It is crucial to have a good phenomenological model of electronic stopping power for modeling 

Q\ t the physics of ion implantation into crystalline silicon. In the spirit of the Brandt-Kitagawa effective 

^\ . charge theory, we develop a model for electronic stopping power for an ion, which can be factorized 

into (i) a globally averaged effective charge taking into account effects of close and distant collisions 
by target electrons with the ion, and (ii) a local charge density dependent electronic stopping power 
for a proton. This phenomenological model is implemented into both molecular dynamics and Monte 
Carlo simulations. There is only one free parameter in the model, namely, the one electron radius 
qq , r° for unbound electrons. By fine tuning this parameter, it is shown that the model can work 

■ successfully for both boron and arsenic implants. We report that the results of the dopant profile 

simulation for both species are in excellent agreement with the experimental profiles measured by 
secondary-ion mass spectrometry (SIMS) over a wide range of energies and with different incident 
directions. We point out that the model has wide applicability, for it captures the correct physics of 
electronic stopping in ion implantation. This model also provides a good physically-based damping 
i' mechanism for molecular dynamics simulations in the electronic stopping power regime, as evidenced 

by the striking agreement of dopant profiles calculated in our molecular dynamics simulations with 
the SIMS data. 
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I. INTRODUCTION 



Ion implantation in semiconductors is an important technology in integrated circuit device fabrication^]. A reliable 
description of as-implanted profiles and the resulting damage is needed for technological development, such as device 
design and modeling, as well as process optimization and control in the fabrication environment. For semiconductor 
devices whose physical dimensions are of order of submicrons or smaller, low implant energies and reduction of thermal 
I/"") ' processing are necessary, resulting in more prominent channeling effects in the as-implanted profiles and less post- 
implant diffusion. At these physical dimensions, it is essential to obtain the two- or three-dimensional details of the 
ever shallower and more compact dopant and damage profiles for_post-implant diffusion simulations. 

Study of the energy loss of channeled particles has a long historyo, for the channeling features can be used to elucidate 
the energy-loss mechanisms. Earlier analytical treatments of the implant profiles based on moment distributions, 
derived from the Lindhard-Scharff-Schiott theory (LSS)S, preclude channeling because of the amorphous nature of 
the targets assumed in the studies. Later, it was realized that, because of the channeling effect, electronic stopping 
. ' power plays a much more significant role in ion implantation into crystalline solids than otherwise would be deduced 
from the applicationQ of the LSS theory to amorphous materials... It is especially true for heavy ion implants at low 
energies, such as arsenic ions in the energy range below 700keV£m. For implantation into silicon, most Monte Carlo 
O | (MC) models are only concerned with boron implants, and have not modeled arsenic implants accurately with an 
electronic stopping power model consistent with that used for boronQtl. As will be shown below, the phenomenological 
model we developed for electronic stopping power can be implemented into a Monte Carlo simulation program for 
both boron and arsenic implants in different channels with equal success over a wide range of implant energies. 

In addition to Monte Carlo simulations with the binary collision approximation (BCA)EI, molecular dynamics (MD) 
incorporating multiple-interactions via many-body potentials can also be used to simulate the behavior of energetic 
ions in amorphous or crystalline silicon. This method is especially applicable at low energies, for which many-body, and 
multiple interactions are increasingly important^! Although it is well known that the BCA is valid for high incident 
energies (~ O.lkeV up to ^MeV, the upper limit is set by relativistic effects), in a cascade, especially initiated by 
a relatively low energy ion, the energy of the ions decrease and eventually reach the lower validity limit of the 
BCA at which many-body effects become impof.tantStj. For crsytals of high symmetry, the BCA can be modified to 



account for simultaneous collisions in channel and MD results can provide good insight into how to successfully 
modify the BCA in this situation. Moreover, MD results can be compared to BCA Monte Carlo simulations and used 
to establish the low energy limits of the binary collision approximation. 
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An extremely important issue in deploying molecular dynamics to model collision processes in covalent and ionic 
solids is how to incorporate energy transfer mechanisms between electrons and ionscj. A good description of dynamical 
processes in energetic collisions, such as initial displacement damage, relaxation processes, and the cooling phase as 
the energy dissipates into the ambient medium, requires a theoretical framework that encompasses all interactions 
between ion-ion, ion-electron, and their interaction with the thermal surroundings. Especially, it should capture 
the noncquilibrium thermodynamic nature of these physical processes involving a wide rangft-ofj-energy scales, from 
a low energy electron-phonon interaction regime to a high energy radiation damage regime .E3~t3. Traditional MD 
simulations can capture the thermal behavior of an insulator. Since they do not take into account coupling between 
the phonons and the conduction electron system, obviously, these simulations underestimate the heat-transfer rate 
for noninsulating materials. In addition to lattice thermal conductivity, the issue of the conductivity due to electrons 
must be addressed. Furthermore, a correct description of the electronic stopping power should be incorporated 
into MD simulations for high energy implantation. For example, in sputtering processes by particle bombardment, 
examination of MD simulations with and without inelastic electronic energy loss has established that, independent 
of the ion's mass or energy, the inelastic electronic energy losses by target atoms within the collision cascade have 
greater influence on the ejected atom yield than the ion's electronic lossesli^l. This is in contrast to the belief that the 
electronic, loss mechanism is important only for cascades initiated by light ions or by heavy ions at high bombardment 
energiestj. Although a convincing experimental verification of the electronic effects in sputtering is still lacking, the 
effects should be relevant to defect production rates, defect mobility and annealing, etcJiSEl Also as shown in Ref. 
pl| , traditional MD simulations produce extremely long channeling tails due to the absence of electronic stopping. In 
order to incorporate the ion-electron interaction into molecular dynamics simulations, a simple scheme was proposed 
by adding a phenomenological term, which describes the inelastic electronic stopping in the high energy radiation 
damage regime, while also capturing the thermal conductivity by coupling low energy ions to a thermal reservoirE3. 
The empirical expression used in Ref. ^2] for the strength of the ion-electron coupling is a function of the local 
electronic density. At the low charge density limit, a density functional result was reproduced^, and at the high 
charge density limit, the linear response results were captured. In the same spirit, we develop a stochastic MD model 
incorporating the electronic stopping power as a damping mechanism. Our model is based on an effective charge 
theoryE-3 with the electronic stopping power factorized into two parts. One is the effective charge of the incident 
ion, which is a globally averaged quantity determined by the average unbound electron density in the medium. The 
other factor is the electronic stopping power for a proton, for which the same local density functional results are used. 
Naturally, our damping mechanism incorporates both regimes, i.e., the electronic stopping regime and the electron- 
phonon interaction regime, into our molecular dynamics simulation, because the inelastic loss for a proton exhibits 
a similar density dependence as prescribed in Ref. ^2], with additional modifications due to the velocity dependence 
of the effective charge. In the present work, however, we emphasize mainly the electron stopping power in the high 
energy regime (~ keV to ^lOOkeV), i.e., the electrons behave as an energy sink. The validity of the model for the 
electronic heat conduction regime will be discussed elsewhere. In the following, for boron and arsenic implants into 
single-crystal silicon in both the channeling and off-axis directions, we will show that a classical MD with the physically- 
based damping mechanism can generate dopant profiles in excellent agreement with experimentally-measured profiles 
obtained by secondary- ion mass spectrometry (SIMS). 

As discussed above, the phenomenological model we have developed for electronic stopping power is successfully 
implemented into both BCA Monte Carlo programs and MD simulations. Wide applicability requires that a model 
be valid for different implant species over a wide range of energies. We emphasize that this electronic stopping model 
is accurate both for boron and arsenic implants, thus providing a crucial test of the generality and validity of the 
model in capturing the correct physics of electronic stopping. 

The paper is organized as follows. In Sec. II, we present the phenomenological model for electronic stopping power 
in detail. Atomic units e = h = m e — I are used throughout the paper unless otherwise specified. In Sec. Ill, 
we briefly discuss differ.eHt electronic stopping models implemented on the versatile BCA Monte Carlo simulation 
platform, MARLOWlilMj-Then the results of the BCA Monte Carlo simulations on a rare-event algorithm enhanced 
UT-MARLOWE platformcj with our electronic stopping model are summarized. In Sec. IV, the results of the MD 
with the inelastic electronic energy loss are presented. In Sec. V, we make closing remarks and point out directions 
for future studies. 



II. THE MODEL 

According to the Brandt-Kitagawa (BK) theoryEi, the electronic stopping power of an ion can be factorized into 
two components based on an effective charge scaling argument. One is the effective charge of the ion (if not fully 
ionized), Z*, which is in general a function of ion velocity v and the charge density of the target p, or equivalently, 
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the one electron radius r s = [3/ (47rp(x))] 1 / 3 ; the other is the electronic stopping power for a proton, S p (v, r s ). In the 
local density approximation, therefore, the total inelastic energy loss AE e of an ion of constant velocity v is 



AE e = J [Zf(«, 



r s )] 2 S p (v,r s )dx, 



(1) 



where the integral is along the ion path. Since the effective charge is a continuous function of electronic density, 
mathematically, it is always possible to find a mean value, r° of r s , such that Eq. (|l|) can be rewritten as 

(2) 



AE e = [ZUv,r°)} 2 J S p (v,r a )dx. 

If the effective charge is a slowly varying function of space, physically, this means that r° describes an average number 
of unbound electrons in the sea and thus can be assumed to determine the Fermi surface. Therefore, we have the 
relation between the Fermi velocity and r°: 
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(3) 



where a = [4 : /(9'k)] 1 ^ 3 . We note that this r° will be the only tunable parameter in our electronic stopping power 
model. 

Next we turn to a simple statistical model for this partially ionized, moving projectile. For an ion with N = Z\ — Q 
bound electrons, where Q is the charge number of the ion of atomic number Z\ , a radially symmetric charge density 
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is used in the BK theory. Here A is the ion size parameter, a function of the fractional ionization, q = (Zi—N)/Zi . The 
total energy of the electrons comes from the sum of the kinetic energy estimated by the local density approximation, 
the electron-electron interaction in the Hartree approximation weighted by a variational parameter A to account for 
correlation, and the Coulomb energy of the electrons in the electric field of the nucleus. A variational approach 
minimizing the total energy leads to the following dependence of the ion size on the ionization fraction q: 



A 



2o (l - q) 2 ' 3 
Zl /3 [l-(l-g)/7] : 



(5) 



where ag = 0.24005. In the BK theory, the generalized Lindhard theory of the electronic stopping in a homogeneous 
electron gas with an electron density n = 3/(47rfg) is used. The total electronic stopping is estimated from the sum 
of the energy loss in soft, distant collisions, i.e., small momentum transfers with target electrons seeing a charge 
qZ\. and the energy loss to the target electrons experiencing increased nuclear interaction in hard, close collisions 
corresponding to large momentum transfers. As extensively discussed in the literature (see, e.g., Ref. HIlll, and 
references therein), it is assumed that the charge state of a proton in a solid is unity. Given an ionization fraction 
q and using the scaling argument for the ratio of ion stopping to the proton stopping at the same velocity, the BK 
theory produces a simple expression for the fractional effective charge of an ionoE3 



7 (f s ) =g + C(f s )(l-g)ln 







where C(r s ) is weakly dependent on the target and has a numerical value of about 1/2. We will set C 
Then, the effective charge is 

Z\ = Z l7 (f s ). 



(6) 



0.5 below. 



(7) 



For our model, using the procedure (Q) outlined above, this dependence of f s is identified with the dependence of the 
mean value r°. Therefore, the effective charge Z* has a nonlocal, i.e., spatially independent, character and depends on 
the Fermi surface. In the above discussion, as can be seen, q is a parameter which is not fixed by the-BK theory. For 
obtaining this ionization fraction, there are velocity and energy criteria originally proposed by BohJH3 and LambEa, 
respectively. Kitagawa also used a statistical argument to justify scaling analyses in terms of the scaling parameter 
v i/(vbZ 2 ^ 3 J^3. Recently, the issue of which stripping criterion can give rise to a better physical understanding has 
been raisedEil. However, in light of the large amount of experimental data employed in Ref. &l\ to extract an ionization 
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scaling consistent with the Brandt-Kitagawa theory, we will use this empirically verified scaling in our model. As 
summarized in Ref. a new criterion in the BK approach is proposeduLJ, i.e., a relative velocity criterion, which 
assumes that the electrons of the ion which have an orbital velocity lower than the relative velocity between the ion 
and the electrons in the medium are stripped off. The relative velocity v r is obtained by averaging over the difference 
between the ion velocity vi and the electron velocity v e under the assumption that the conduction electrons are a free 
electron gas in the ground state, therefore, whose velocity distribution is isotropic. Performing a further averaging of 
v e over the Fermi sphere leads tcEa 



v r 
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for v\ > vf, 
for v\ < vp- 



(8) 
(9) 



For the ionization scaling, a form of the Northcliffe typeEa is then assumed for the scaling variable, i.e., the reduced 
relative velocity: 
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where vb is the Bohr velocity and vb 
used in Ref. [27| to determine 



1 in our units. The extensive experimental data for ions 3 < Z\ < 92 are 



q = 1 - exp[-0.95(y r - 0.07)]. 



(11) 



ital data along the fit is presented. 
1 The accuracy level of Eq. ( [lT] ) is 



In Ref. [27| an ionization scaling fit with even tighter bunching of the experime 
However, this approach entails a much more involved computational proceduret 

adequate for our present purposes. .— ■ 
In our model, the electronic stopping power for a proton is derived from a nonlinear density-functional formalisms. 
In the linear response theory, the energy loss per unit path length of a proton moving at velocity v in the electron 
gas is obtained by Ritchieo 
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(12) 



using an approximation to the full random-phase approximation dielectric function, which amounts to the exponential 
screening potential around the ion induced by density fluctuations of the electrons .-The nonlinear, density- functional 
calculation based on the formalism of Hohenberg and Kohn, and Kohn and ShamEJO has been performedE§E30 to 
obtain the charge density and scattering phase shifts for the conduction band as a function of energy self-consistently. 
The final stopping power for a proton is obtained via 
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(Z + l)sin 2 [6i(E F ) - 6i +1 {E F )] , 



(13) 



where Si(Ef) is the phase shift at the Fermi energy for scattering of an electron of angular momentum I and fcf is 
the Fermi momentum!^. As shown in Ref. p7| , ^9p0| , comparison with expermental data cLemonstrates that the density 
functional treatment provides an improvement over the linear response (dielectric) resultEJ, which underestimates the 
stopping powers. In our implementation, the result of the nonlinear, density functional formalism for the electronic 
stopping power for a proton is used, which can be expressed as 



S '""'-' J ( ^f) (V(r - ) ' 



(14) 



where, for computational convenience, the correction factor G(r s ) takes the form 

G(r s ) = 1.00 + 0.717r s - 0.125r^ - 0.0124?^+ 0.00212r^ (15) 

for r s < 6. We note that a different correction factor was used in Refs. p2p^ , which does not have the following 
desired behavior for r, < 1. Since the density functional result converges to the Ritchie formula as r s decreases 
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towards values sufficiently small compared to unitya, this requires that the correction factor smoothly tend to unity 
as r s — > 0. Obviously, the above G(r s ) possesses the correct convergence property. 

The last ingredient needed for our model is the charge_d,istribution /o(x) for silicon atoms in the crystal. We 
use the solid-state Hartrcc-Fock atomic charge distributionEZl, which is spherically symmetric due to the muffin-tin 
construction. In this approximation, there is about one electron charge unit (0.798 electrons for Si) left outside the 
muffin-tin. This small amount of charge can be either distributed in the volume between the spherical atoms, resulting 
in an interstitial background charge density 0.119e/A , or distributed between the maximal collision distance used in 
Monte Carlo simulations and the muffin-tin radius (see details below). 



III. BCA MONTE CARLO SIMULATION RESULTS 



First, in comparison with other electronic stopping models used in Monte Carlo simulations based on the MAR- 
LOWE platform, we stress that in our model the effective charge is a nonlocal quantity, neither explicitly dependent 
on the impact parameter nor on the charge distribution, and the stopping power for a proton depends on the-Lpcal 
charge density of the solid. A purely nonlocal version of the BK theory was implemented into MARLOWEej, in 
which both the effective charge and the stopping power for a proton depend on a single nonlocal parameter, namely, 
the averaged one electron radius. Its results demonstrated that energy loss for well-channeled ions in the keV region 
has high sensitivity to the one-electron radius in the channel. Ii,.wp,s pointed out that a correct density distribution 
is needed to account for the electronic stopping in the channefiea. Later, ^qpirely local version of the BK theory 
was developed to take into account the charge distribution of the electronsMr 4 !. Comparison with other electronic 
stopping models, such as Lindhard and Scharffl, Firso'vCj, ag d th e above nonlocal implementations, showed a marked 
improvement in modeling electronic stopping in the channeflSci Good agreement between simulated dopant profiles 
and the SIMS profiles for boron implants into (100) single crystal silicon was obtained. However, this purely local 
implementation of the BK theory did not successfully model the electronic stopping for the boron implants into the 
(110) axial channel and arsenic implants-as noted in Ref. 0. 

In the present work, UT-MARLOWEca was selected as the platform for our electronic stopping model implemen- 
tation. UT-MARLOWE is an extension of the MARLOWE code for simulating the behavior of energetic ions in 
crystalline materialsBlMl It has bejen enhanced with: (i) atomic pair-specific interatomic potentials for B-Si, B-O, 
As-Si, As-O for nuclear stoppingsEZI, (ii) variance reduction algorithm implemented for rare events, (iii) important 
implant parameters accounted for, e.g., tilt and rotation angles, the thickness of the native oxide layers, beam diver- 
gence, and wafer temperature, etc. In our simulations, we have turned off certain options, such as the cumulative 
damage model in the UT-MARLOWE code, which is a phenomenological model to estimate defect production and 
recombination rates. Individual ion trajectories were simulated under the BCA and the overlapping of the damage 
caused by different individual cascades was neglected. In order to test the electronic stopping model we also used low 
dose (10 13 /cm 2 ) implants so that cummulative damage effects do not significantly complicate dopant profiles!! Also, 
for the simulation results we report below, 16A native oxide surface layer, 300K wafer temperature were used. Th-p 
maximum distance for searching a collision partner is 0.35 lattice constant, the default value in the UT-MARLOWEH. 
The excess charge outside the muffin-tins is distributed in the space between this maximum collision distance and the 
muffin-tin radius. In the simulation, the electronic stopping power is evaluated continuously along the path the ion 
traverses through regions of varying charge density, i.e., the energy loss is given by 

AE e = f [Zi7(«i,0] 2 S p («i,r a (x))da;. (16) 

J ion path 

In the simulations, the free parameter r° was adjusted to yield the best results in overall comparison with the 
experimental data. The value r° — 1.109A was used for both boron and arsenic ions for all energies and incident 
directions. This value is physically reasonable for silicon. Note that the unbound electronic density in silicon with 
only valence electrons taken into account will give rise to a value of 1.06lA for r s . The fact that our r° value is greater 
than 1.06lA indicates that not all valence electrons participate in stopping the ion as unbound electrons. 

We display the Monte Carlo dopant profile simulation results as follows. We note in passing that the lower and 
upper limits of energy used in our simulations are determined by the energy range of the SIMS data available to us. 

In Figs, [lj U and||, we show boron dopant profiles for the energies 15keV, 35keV, and 80 keV along (100), (110), and 
the off-axis direction with tilt = 7° and rotation = 30°, respectively. It can be seen that the overall agreement with 
the SIMS data is excellent. In Fig. |IJ the simulations show a good fit for the cutoff range. In the high energy regime, 
the simulated distribution shows a slightly peaked structure. This can be attributed to a strong channeling due to 
insufficient scatterings in the implants. We have noticed that by increasing, e.g., the native oxide layer thickness, the 
peak can be reduced. For the (110) channeling case, the distribution indicates a possibility that the total electronic 
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stopping power along the channel is a little too strong at the high energy end. However, it should be kept in mind that 
for this channel, the UT-MARLOWE model becomes sensitive to the multiple collision parameter which is employed 
as an approximate numerical correction for the effect of multiple overlapping nuclear encounters. It is not clear how 
to separate the contributions from these two different sources. ■— . 

For comparison, in Fig. |], we also display a low energy (5keV) implant casecj. Again the agreement for both 
the channeling and off-axis directions is striking (thin lines without symbols). In order to illustrate the importance 
of electronic stopping power at this low energy for boron implants, we have used an artificially reduced electronic 
stopping power, i.e, multiplying AE e in Eq. ( |l6| ) by a factor of 1/10 in the simulation, to genearte dopant profiles in 
the channeling and off-axis directions. From Fig. ^, evidently, it can be concluded that, for boron implants even in this 
low energy regime, electronic stopping power has significant influence on the channeled tail of the dopant distribution 
and on the cutoff range for both channeling and off-axis directions. 

Figs. H and || show arsenic dopant profiles for energies ranging from 15keV to 180 keV along (100), and the off-axis 
direction with tilt = 8° and rotation = 30° , respectively. It can readily be concluded that our electronic model works 
successfully with arsenic as well as boron implants into crystalline silicon. For comparison, a case of arsenic implant 
into amorphous silicon is also shown in Fig. [7]. The implant energy is 180kcV. The effect of electronic stopping, which 
shows clearly in the long sloping channeling tail in the crystalline counterpart (see Fig. ^|), is less prominent for the 
amorphous case (Fig. [?]). 

To examine the role that electronic stopping power has on arsenic implants in the low energy regime (~ 5keV), we 
again simulated arsenic dopant profiles with the artificially reduced electronic stopping power. For these low energy 
implants, the oxide layer thickness 3A were used in the BCA simulations on account of the fact that the wafers used 
for these implants were treated by dilute HF etch for 30 seconds, then implanted within 2 hours to prevent native oxide 
regrowthcll. In Fig. ^|, we show that our electronic stopping power model is successful in both the (100) channeling 
and off- axis directions (thin lines without symbols). Clearly, the artificial reduction of electronic stopping power 
leads to incorrect dopant distributions and cutoff ranges for both the channeling and off-axis directions, although the 
deviations indicate a less significant contribution from electronic stopping for arsenic implants than for boron implants 
at the energy 5keV. However, the deviation in the cutoff range due to the electronic stopping power reduction for the 
channeling case is still significant. Obviously, this reinforces the conclusion that, for channeling implants even at low 
energies, electronic stopping is not negligible. 

In summary, the above results demonstrate clearly that our electronic stopping power indeed captures the correct 
physics of the electronic stopping for ion implants into silicon over a wide range of implant energies. 

IV. MD SIMULATION RESULTS 

We have also used classical molecular dynamics simulation to study the electronic stopping power as one of the 
damping mechanisms in the high energy regime, as discussed above. Here we demonstrate that experimental data, 
such as SIMS, can be used to test the validity of this_physically-based damping model. The interaction between silicon 
atoms are modeled by Tersoff 's empirical potentialEj: 

E = \ E H r v) \y^ r a) - b Mnj)} , (17) 

where f(rij) is a cutoff function that restricts interactions to nearest neighbors, Vfj(ry) and Vyj^ry) are pair terms, and 
bij is a many-body function that can be regarded as an effective Pauling bond order. _We have modified the repulsive 
part of the Tersoff potential by splinning to the ZBL universal potential at close-range .ED. The ZBL universal potential 
is also used to model the ion-silicon interactions. In our full MD simulations for the low dose implantation, the lattice 
temperature was initialized to 300K and the above electronic stopping model was applied to all the atoms. The only 
modification required for implementation in MD is to take into account the contributions from multiple silicon atoms 
to the local electron density, while ensuring that the background electron density is only counted once. For each 
individual cascade, all recoils and the accumulation of damage in the ion path are taken into account. Using the 
parameter value 1.109 for r° from the comparison of BCA Monte Carlo simulation results with the SIMS data, we 
have simulated the implantation of low energy boron and arsenic ions into the Si {100}(2 x 1) surface at energies 
between 0.5keV and 5keV, with both channeling and off-axis directions of incidence. We mention here that, for the 
(100) channeling case up to 0.16 keV (32%) of 0.5 keV boron implant energy and 0.64 keV (13%) of 5keV arsenic 
implant energy are lost via electronic stopping in our simulations. Simulations were terminated when the total 
energy of the ion became less than 5eV, giving typical simulation times of around 0.2ps. Figs. [)[ [lO] and [ll] show 
the calculated dopant concentration profile for various energies and directions. Each MD profile is generated by a 
set of between 500 and 1300 individual ion trajectories. Also shown are the profiles obtained using the modified 
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UT-MARLOWE BCA code described in Sec. III. Obviously, the MD calculation results are in very good agreement 
with the experimental data, and with the BCA results. This demonstrates that our electronic stopping power model 
provides a good physically-based damping mechanism for MD simulations of ion implantation. 



V. CONCLUSION 



We have developed a phcnomcnological electronic stopping power model for the physics of ion implantation. It has 
been implemented into MD and BCA Monte Carlo simulations. SIMS data have been used to verify this model in 
the MD and BCA Monte Carlo platforms. This model has only one free parameter, namely, the one electron radius 
of unbound electrons in the medium. We have fine tuned this parameter to obtain excellent results of dopant profiles 
compared with SIMS data in both MD and BCA Monte Carlo simulations. We emphasize that this model with a 
single parameter can equally successfully model both boron and arsenic implants into silicon over a wide range of 
energies and in different channeling and off-axis directions of incidence. This versatility indicates wide applicability of 
the model in studies of other physical processes involving electronic stopping. As a more stringent test of the model, 
it should also be applied to implantation of species other than boron and arsenic. Using arsenic implantation as an 
example, we have also addressed the issue of how significant electronic stopping is for heavy ions in a low energy 
regime. For instance, to achieve a good quantitative understanding, we still have to take into account the physics of 
electronic stopping for arsenic implants at 5keV. 

As discussed above, it is important to incorporate ion-electron couplings into MD simulations in both the high energy 
radiation damage regime and the low energy electron-phonon interaction regime. We have demonstrated that this 
model provides a crucial piece of physics in MD simulations for modeling energetic collisions in the electronic stopping 
power regime. The agreement of the simulated dopant profiles with the SIMS data shows that the incorporation 
of this physically-based damping term into MD simulations is a phenomenologically reliable approach in the regime 
concerned. Under way is an investigation of whether it can be used as a good phenomenological model for electron- 
phonon coupling in the low energy regime. This agreement also suggests that MD can be used to generate dopant 
profiles for testing against the low energy BCA results when experimental data is not available. Furthermore, MD 
simulations incorporating this physically-based damping mechanism can provide valuable insight into how to modify 
the binary collision approximation. This will enable the validity of the Monte Carlo simulation to be extended 
further into the lower energy regime, while not destroying computational efficiency required in realistic simulation 
environments. 
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VII. FIGURES 
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FIG. 1. Boron implant profiles for the (100) direction with energies ranging from 15keV to 80keV. Zero tilt and rotation 
angles. The thick lines are SIMS data. 
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FIG. 2. Boron implant profiles for the (110) direction with energies ranging from 15keV to 80keV. Zero tilt and rotation 
angles. The thick lines are SIMS data. 
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FIG. 3. Boron implant profiles for the (100) direction with the tilt = 7° and rotation = 30°. The thick lines are SIMS data. 
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FIG. 4. Boron implant profiles for the (100) channeling and for the off-axis direction with the tilt = 7° and rotation = 7°. 
The implant energy is 5keV. Thick lines: SIMS data; Thin lines without symbols: BCA Monte Carlo simulation with full 
electronic stopping power; Circles: BCA Monte Carlo simulation with the artificially reduced electronic stopping power for the 
channeling case; Triangles: the corresponding case for the off-axis direction (see text). 
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FIG. 5. Arsenic implant profiles for the (100) direction with energies ranging from 15keV to 180keV. Zero tilt and rotation 
angles. The thick lines are SIMS data. 
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FIG. 6. Arsenic implant profiles for the (100) direction with the tilt = 8° and rotation = 30°. The thick lines are SIMS data. 
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FIG. 7. Arsenic implant profile into amorphous silicon with the implant energy being 180keV. The thick line is SIMS data. 



E 



10 ,s 



10 ,s 























\ \ 


<100> channel, angel=0,0 








V '■ 






angle=10°,22° \ V 


i'UH I. I V\. ,,>>, 



U , i , i S * L-i 1— a Li I \ . . > , ^ , | 

0.0 0.2 0.4 0.6 0.8 1.0 

DEPTH (10* ANGSTROM) 



FIG. 8. Arsenic implant profiles for the (100) channeling and for the off-axis direction with the tilt = 10° and rotation 
= 22°. The implant energy is 5keV. Thick lines: SIMS data; Thin lines without symbols: BCA Monte Carlo simulation with 
full electronic stopping power; Circles: BCA Monte Carlo simulation with the artificially reduced electronic stopping power for 
the channeling case; Triangles: the corresponding case for the off-axis direction (see text). 
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FIG. 9. Molecular dynamics simulation: boron implant profiles into the (100) channel with tilt = 10°, rotation = 22°. 
Implant energy is 5keV. Thick line: SIMS data. Thin line: MD simulation. Dotted line: BCA Monte Carlo simulation. 
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FIG. 10. Molecular dynamics simulation: arsenic implant profiles into the (100) channel with zero tilt, rotation angles. 
Implant energy is 5keV. Thick line: SIMS data. Thin line: MD simulation. Dotted line: BCA Monte Carlo simulation. 
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FIG. 11. Molecular dynamics simulation: arsenic implant profiles into the (100) channel with tilt = 10°, rotation 
Implant energy is 5keV. Thick line: SIMS data. Thin line: MD simulation. Dotted line: BCA Monte Carlo simulation. 
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